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1. Introduction 



The semiclassical analysis has proven to be a very useful tool in the field of "quantum 
chaos" as well as in many other fields. Nevertheless, one should bear in mind, that 
it only approximates the true quantal solution. Thus, it is imperative to know what 
are the errors which are inherent to the semiclassical approximation, and whether they 
could be considered as sufficiently small for the problem at hand. 

We shall focus our attention on one particular application of the semiclassical 
approximation: The calculation of the energy spectra of classically chaotic systems. 
The analytical tool that is used for this purpose is the semiclassical Gutzwiller trace 
formula H which expresses the quantum spectral density in terms of classical quantities, 
and in particular the actions and stabilities of classical periodic orbits. The trace 
formula was used, among other things, to explain and discuss spectral statistics and their 
relation to the universal predictions of Random Matrix Theory (RMT) ||, [|. However, 
a prerequisite for the use of the semiclassical approximation to compute short-range 
statistics is that it is able to reproduce the exact spectrum within an error which is 
comparable to or less than the mean level spacing! This is a demanding requirement, 
and quite often, the ability of the semiclassical approximation to reproduce precise levels 
for high-dimensional systems is doubted, and on the following grounds. The mean level 
spacing depends on the dimensionality (number of freedoms) of the system, and it is 
C?(ft d )@. Gutzwiller uses an argument by Pauli [||] to show that in general the 
error margin for the semiclassical approximation scales as 0(h 2 ) independently of the 
dimensionality. Applied to the trace formula, the expected error in units of the mean 
spacing, which is the figure of merit in the present context, is therefore expected to 
be 0{h 2 ~ d ). We shall refer to this as the "traditional estimate". It sets d = 2 as a 
critical dimension for the applicability of the semiclassical trace formula and hence for 
the validity of the conclusions which are drawn from it. The few systems in d > 2 
dimensions which were numerically investigated display spectral statistics which adhere 
to the predictions of RMT as accurately as their counterparts in d = 2 || |7|, || . Thus, 
the "traditional estimate" cannot be entirely correct in the present context, and we shall 
explain the reasons why it is inadequate when we discuss other error bounds in the next 
section. 

It is rather surprising that the problem of the accuracy of the semiclassical trace 
formula was rather rarely touched upon in the literature. Gutzwiller quotes the 
"traditional" estimate of 0(h 2 ~ d ) which was discussed above |lj, ||. Gaspard and Alonso 
|TIJ and Alonso and Gaspard [[□]] derived explicit (generic) h corrections for the periodic 
orbit terms in the trace formula, but did not investigate their effect on the semiclassical 
accuracy of energy levels. Diffraction corrections were discussed in the context of the 
trace formula by Vattay, Wirzba and Rosenqvist and by Primack et al. [I3|. Also 
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in these works the focus is on the corrections to individual periodic orbit terms rather 



than on the overall effect on energy levels. Boasman |TJ| studied the accuracy of the 



Boundary Integral Method (BIM) ]T5[ for 2D billiards in the case that the exact kernel is 
replaced by its semiclassical asymptotic approximation. He finds that the resulting error 
is of the same magnitude as the mean spacing, which is consistent with the traditional 
estimate. However, Boasman's work does not refer directly to the trace formula and 
to the periodic orbits contributions. The works of Bleher ]n| |T7| and of Prosen and 



Robnik |T8| discuss the accuracy of the semiclassical approximation in the integrable 
case. 

The purpose of this work is to address the subject of the accuracy of the semiclassical 
trace formula conceptually, theoretically and numerically. We shall be interested in 
particular in the dependence of the semiclassical error on the dimension of the system. 
To do so, we shall have to start by developing the basic concepts and define the measures 
we use for a quantitative estimate of the spectral error (section §). The accuracy of the 
semiclassical approximation of the quantal energy spectrum will be then studied via 
the dual classical spectrum of actions and stabilities of periodic orbits (time spectrum). 
This will enable us to use our data base of quantum levels and periodic orbits for the 
Sinai billiards in two and three dimensions for a direct evaluation of the semiclassical 
error (section |3|). We shall summarize the paper and discuss a few relevant points in 
section £|. 

2. Measures of the semiclassical error 

In order to define a proper error measure for the semiclassical approximation of the 
energy spectrum one has to clarify a few issues. In contrast with the EBK quantization 
which gives an explicit formula for the spectrum, the semiclassical spectrum for chaotic 
systems is implicit in the trace formula, or in the semiclassical expression for the spectral 
determinant. To extract the semiclassical spectrum we recall that the exact spectrum, 
{E n }, can be obtained from the exact counting function: 

oo 

N(E) = £ G(E - E n ) , (1) 



n=l 

by solving the equation 

N(E n )=n-±, n=l,2,... . (2) 

In the last equation, an arbitrarily small amount of smoothing must be applied to the 
Heavyside function. In complete analogy, one obtains the semiclassical spectrum {Eff} 
as 

N sc (E s n c )=n-^, n = l,2,... , (3) 
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where N sc is a semiclassical approximation of N. Note that iV sc with which we start 
is not necessarily a sharp counting function. However, once {E^} is known, we can 
"rectify" the smooth N sc into the sharp counting function iV# ||: 

oo 

N*(E) = £ &{E - K) . (4) 

n=l 

The simplest choice for N sc is the Gutzwiller trace formula truncated at the 
Heisenberg time, which is what we shall use in the present paper. Alternatively, one 
can start from the regularized Berry-Keating Zeta function ( S c{E) ||20|| , and define 
^sc = (V 71 ") I m 1°8 Csc{E + iO), in which case N sc = iV#. 

Next, in order to define a quantitative measure of the semiclassical error, one should 
establish a correspondence between the quantal and the semiclassical levels, namely, one 
should identify the semiclassical counterparts of the exact quantum levels. In classically 
chaotic systems, for which the Gutzwiller trace formula is applicable, the only constant 
of the motion is the energy. This is translated into a single "good" quantum number in 
the quantum spectrum, which is the ordinal number of the levels when ordered by their 
magnitude. Thus, the only correspondence which can be established between the exact 
spectrum {E n } and its semiclassical approximation, {E^}, is 

En — K . (5) 

This is to be contrasted with integrable systems, where it is appropriate to compare the 
exact and approximate levels which have the same quantum numbers. 

The scale on which the accuracy of semiclassical energy levels should be measured 
depends, in general, on the problem at hand. The most natural choice, however, is the 
mean level spacing (d(E))~ l where d is the smooth density of states. The semiclassical 
error of the n'th level is therefore measured by |14| : 

e n = d(E n ) (E n - E^) . (6) 

A more useful and significant measure is the average of e n over an energy interval AE 
centered at E, which contains a large number of quantum energies, but which must 
be so small that both the classical dynamics and the mean density of states remain 
approximately constant. This energy averaging will be denoted by triangular brackets 
(•) e in the sequel. If the semiclassical mean density d sc agrees with J to a high precision, 
then obviously e(E) = (e n )E = 0. In this case e(E) is not a useful error measure. For 
billiard systems this is always the case, since the mean spectral density can be written as 
an asymptotic series with explicitly known coefficients [[21], ^2], |23|. For general systems, 
only the leading Weyl term is explicitly known. Two appropriate measures which are 
sensitive to the accuracy of the fluctuating parts of the level densities are the mean 
absolute error: 

e ( - 1 \E) = {d(E n )\E n -E s n c \) E (7) 
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and the variance: 

e (2) (£) = < (d(E n )(E n -E-)f) E . (8) 

Having defined the spectral error measures, let us apply them and try to get some 
estimates of the semiclassical error. In the introductory section we have mentioned 
the "traditional" estimate of the semiclassical error. Gutzwiller |9] shows, based on 
Pauli |J , that the inherent error in the semiclassical ( Van- Vleck) approximation of the 
quantal time propagator scales like 0(h 2 ). Since the energy levels are the temporal 
Fourier components of the propagator, it is plausible to assume that they have the same 
degree of accuracy: 

E n - E s n c = 0(h 2 ) "traditional". (9) 

(Strictly speaking, this is an upper bound.) The mean density of energy levels, for a 
general d-dimensional system is asymptotically given by Weyl's formula [|]]: 

d( E ) « = °^ d ) ( 10 ) 

where u(E) is the measure of the energy surface in the classical phase space. Hence, 

traditional = G -2-dx f COnst , d = 2 ag ft — > 0. (11) 

[oo , a > 3 

That is, the semiclassical approximation is (marginally) accurate in 2 dimensions, but 
it fails to predict accurate energy levels for 3 dimensions or more. 

One may get a different estimate of the semiclassical error, if the Gutzwiller Trace 
Formula (GTF) is used as a starting point. Suppose that we have calculated N sc to a 
certain degree of precision, and we compute from it the semiclassical energies using (§). 
The quality of this approximation can be estimated if the leading corrections AiV sc are 
also included and the resulting energy differences 5 n are evaluated. We thus need to 
consider: 

N sc (E* n c + S n ) + AN sc (E s n c + 8 n ) = n-~. (12) 

Combining (|3p and ( |T^ ) we get to first order in S n : 

^ AN SC (E-) ^ AN SC (E-) 
n ~ dN sc (E%)/dE ~ d{E%) • 1 > 

In the above we assumed that the fluctuations of iV sc around its average are not very 
large. Thus, 

e GTF « d(E%)5 n « AiV sc (^ c ). (14) 

Let us apply the above formula and consider the case in which we take for iV sc its 
mean part N, and that we include in N terms of order up to (and including) h~ m , m < d. 
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For AiV sc we use both the leading correction to N and the leading order periodic orbit 
sum which is (formally) of order h°. Hence, 

4 TF = 0(h- m+l ) + O(h ) = O (h^H-m+l, 0) j _ ( 15 ) 

We conclude, that approximating the energies only by the mean counting function N 
up to (and not including) the constant term, is already sufficient to obtain semiclassical 
energies which are accurate to O(h ) = 0(1) with respect to the mean density of states. 
Note again, that no periodic orbits were included in N sc . Including less terms in N will 
lead to a diverging semiclassical error, while more terms will be masked by the periodic 
orbit (oscillatory) term. One can do even better if one includes in N sc the smooth terms 
up to and including the constant term (O(h )) together with the leading order periodic 
orbit sum which is formally also O(h ). The semiclassical error is then: 



e — = Q{h l ). (16) 



,GTF 
-po 



That is, the semiclassical energies measured in units of the mean level spacing 
are asymptotically accurate independently of the dimension! This estimate grossly 
contradicts the "traditional" estimate (|TT|) and calls for an explanation. 

The first point that should be noted is that the order of magnitude (power of h) of 
the periodic orbit sum, which we considered above to be O(h ), is only a formal one. 
Indeed, each term which is due to a single periodic orbit is of order O(% ). However the 
periodic orbit sum absolutely diverges, and at best it is only conditionally convergent. To 
give it a numerical meaning, the periodic orbit sum must therefore be regularized. This 
is effectively achieved by truncating the trace formula or the corresponding spectral £ 
function [^4], ^5], |20], |26j. However, the truncation cutoff itself depends on h. One can 
conclude, that the simple minded estimate ([II]) given above is at best a lower bound, 
and the error introduced by the periodic orbit sum must be re-evaluated with more 
care. This point will be dealt with in great detail in the sequel, and we shall eventually 
develop a meaningful framework for evaluating the magnitude of the periodic orbit sum. 

The connection and the disparity between the "traditional" estimate of the 
semiclassical error and the one based on the trace formula can be further illustrated 
by the following argument. The periodic orbit formula is derived from the semiclassical 
propagator K sc using further approximations [|]]. One thus wonders, how can it be that 
further approximations of K sc actually reduce the semiclassical error from (|TTD to fll6|)? 
The puzzle is resolved if we recall, that in order to obtain e^ F above we separated 
the density of states into a smooth and an oscillating parts, and we required that the 
smooth part is accurate enough. To achieve this, we have to go beyond the leading 
Weyl's term and to use specialized methods to calculate the smooth density of states 
beyond the leading order. These methods are mostly developed for billiards [El, 22, 23[ . 
In any case, to obtain e^ F we have added an additional information which goes beyond 
the leading semiclassical approximation. 
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A straightforward check of the accuracy of the semiclassical spectrum using the error 
measures e,e^\e^ 2 ' is exceedingly difficult due to the large number of periodic orbits 
needed because of the exponential proliferation in chaotic systems. The few cases where 
such tests were carried out involve 2D systems and it was possible to check only the 
lowest (less than a hundred) levels (e.g. [E7|, The good agreement between the exact 
and the semiclassical values confirmed the expectation that in 2D the semiclassical error 
is small. In 3D, the topological entropy is typically much larger [T9|, [7]], and the direct 
test of the semiclassical spectrum becomes prohibitive. 

Facing with this grim reality, we have to introduce alternative error measures which 
yield the desired information, but which are more appropriate for a practical calculation. 
We construct the measure: 

5®(E) = (jN(E)-N#(E)\ 2 ) E . (17) 

As before, the triangular brackets indicate averaging over an energy interval AE about 
E. We shall now show that 8^ 2 ' faithfully reflects the deviations between the spectra, 
and is closely related to and . Note, that the following arguments are purely 
statistical and apply to every pair of staircase functions. 

Suppose first, that all the differences E^ — E n are smaller than the mean spacing. 
Then, \N — A#| is either or 1 (see figure HJ) and hence \N — iV#| = \N — A#| 2 . 
Consequently, 

5 {2) (E) » (\N(E) - N*(E)\) e , small deviations. (18) 

However, the right hand side of the above equation (the fraction of non-zero 
contributions) equals e^. Thus, 

« , small deviations . (19) 

If, on the other hand, deviations are much larger than one mean spacing, the typical 
horizontal distance d\E — E n \ should be comparable to the vertical distance \N — N#\, 
and hence, in this limit 

#(2) ~ e (2) ? i arge deviations . (20) 

Therefore, we expect to interpolate between e'- 1 -' and e*- 2 -* throughout the entire 
range of deviations. This behavior is indeed observed in a numerical test which was 
performed to check the above expectations. We considered the unfolded exact spectrum 
(normalized to unity mean spacing [pSfl ) of the 3D Sinai billiard {X n } and created from 
it a synthetic spectrum by adding a random variable with mean and variance a 2 : 

X° = X n + X ran dom(0, a). (21) 

The {A^"} has also a unity mean spacing and it is meant to imitate a semiclassical 
spectrum. After sorting {X°} we calculated the measures e> x \e?> and 5^ as functions 
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Figure 1. Illustration of \N(E) - N*(E)\ for small deviat ions between quantum and 
semiclassical energies: <C dT 1 = A. The quantum staircase N(E) is denoted by 
the full line and the semiclassical staircase N#(E) is denoted by the dashed line. The 
difference is shaded. 



of a. The results are shown in figure [|, and they verify the estimates (fL9|, |20|) in 
the appropriate limits. The numerical test reported in figure |] demonstrates another 
attractive feature of the measure 5^: It is completely equivalent to e^> when the 
spectral counting functions are replaced by their smooth counterparts, provided that 
the smoothing width is of the order of 1 mean level spacing and the same smoothing is 
applied to both counting functions. That is, 

"smooth ~ fc 

for all deviations. (In fact, for small deviations there is a proportionality factor, but 
it can be set to 1 if an appropriate smoothing is used.) In testing the semiclassical 
accuracy, this kind of smoothing is essential and will be introduced by truncating the 
trace formula at the Heisenberg time £h = hd. These properties of the measure 
and its complete equivalence to for smooth counting functions, renders it a most 
appropriate measure of the semiclassical error. 

We now turn to the practical evaluation of S^ 2 \ To affect the energy averaging, 
we choose a positive window function w(E' — E) which has a width AE near E and 
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Figure 2. The spectral measures e 1 ), e*- 2 -*, S 1 - 2 ^ and and their asymptotics for a wide 
range of spectral deviations a. The random distribution was Gaussian, and similar 
results were obtained also for a uniform distribution. Note the logarithmic scale. 



is normalized by dE'w 2 (E') = 1. It falls off sufficiently rapidly so that all the 
expressions which follow are well behaved. Construct the following counting functions 
that have an effective support on an interval of size AE about E: 



N(E';E) =w(E'-E)N(E') 
N*{E'-E)=w{E>-E)N*{E>) 



(23) 
(24) 



At this stage N and iV* are still sharp staircases, and we note that the multiplication 
with w preserves the sharpness of the stairs (it is not a convolution!). We now explicitly 
construct 5^(E) as: 



5 (2 \E) = / dE' N(E'; E) - N#(E'; E) 

J — oo 



-oo 

-oo 



/~r°° 2 
dE' N{E')-N*{E') w 2 {E'-E) 



(25) 



To construct <^ ooth we need to smooth N, N# over a scale of order of one mean spacing. 
This can be done e.g. by replacing the sharp stairs by error functions. As for 7V#, we 
prefer to simply replace it with the original N sc , which we assume to be smooth over one 
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a 



Figure 3. Comparison of the measures <$sJ ooth and for a wide range of spectral 
deviations a and for Gaussian distribution. The difference is also plotted since the 
curves almost overlap. Note the logarithmic scale. 



mean spacing. That is, we suppose that N sc contains periodic orbits up to Heisenberg 
time. Hence, 

€!ooth(£) = / dtf iV sm °° th (E>) - N SC (E')\ w\E> - E). (26) 

J — oo 



A comment is in order here. Strictly speaking, to satisfy ( p2| ) we need to apply the same 
smoothing to N and to iV#, and in general jy*' 811100 * 11 ^ N sc , but there are differences 
of order 1 between the two functions. However, since our goal is to determine whether 
the semiclassical error remains finite or diverges in the semiclassical limit h —>■ 0, we 
disregard such inaccuracies of order 1. If more accurate error measure is needed, then 
more care should be practiced in this and in the following steps. 
Applying Parseval's theorem to ( |26| ) we get: 

<5 s ( 2ooth(£) = jrl di D(t; E) - D sc {t; E) (27) 



oo 



where 



1 r+oo 

D(t;E) = -y= I dE' N smooth (E';E)exp(iE't/h) (28) 

v27T J-oo 
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D sc (t; E) = -= dE'N sc (E'; E) exp(iE't/h). (29) 

V Z7T J—oo 

We shall refer to D, D sc as the (regularized) quantal and semiclassical time spectra, 
respectively. This name can be justified by invoking the Gutzwiller trace formula and 
expressing the semiclassical counting function as a mean part plus a sum over periodic 
orbits. We have: 

N SC (E) = N(E) + £ smiSjiEyh - vj-k/2] , (30) 



where Aj = Tj/(7rhrj^J\I — Mj\) is the semiclassical amplitude of the j'th periodic orbit, 
and Tj, Sj,Uj, Mj,rj are its period, action, Maslov index, monodromy and repetition 
index, respectively. Then, the corresponding time spectrum reads: 

D sc (t;E)paD(t;E) (31) 

+ ^i^ej {^ n)[Et+Sj{E)] ^ + urnm- 

In the above, the Fourier transform of w is denoted by w. It is a localized function 
of t whose width is At pa h/AE. The sum over the periodic orbits in D sc therefore 
produces sharp peaks centered at times that correspond to the periods Tj, hence the 
name "time spectrum" . The term D corresponds to the smooth part and thus is sharply 
peaked near t = 0. To obtain ( |3~T{ ) we expanded the actions near E to first order: 
Sj(E') pa Sj(E) + (E' — E)Tj(E). We note in passing, that this approximate expansion 
of Sj can be avoided altogether if one performs the Fourier transform over h~ x rather 
than over the energy. This way, an action spectrum will emerge, but also here the action 
resolution will be finite, because the range of hT x should be limited to the range where 
d(E; h) is approximately constant. It turns out therefore, that the two approaches are 
essentially equivalent, and for billiards they are identical. 

The manipulations done thus far were purely formal, and did not manifestly 
circumvent the difficult task of evaluating 5 sr n 00th . However, the introduction of the 
time spectra and the formula (p7| ) put us in a better position than with the original 
expression ( p5|) . The advantages of using the time spectra in the present context are 
the following: 

• The semiclassical time spectrum D sc (t;E) is absolutely convergent for all times 
(as long as the window function w is well behaved, e.g. it is a Gaussian). This 
statement is correct even if the sum (31) extends over the entire set of periodic 
orbits! This is in contrast with the trace formula expression for iV sc (and therefore 
N sc ) which is absolutely divergent if all of the periodic orbits are included. 
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• Time scale separation: As we noted above, the time spectrum is peaked at times 
that correspond to periods of the classical periodic orbits. This allows us to 

(2) 

distinguish between various qualitatively different types of contributions to Smooth- 
We shall now pursue the separation of the time scales in detail. We first note, that due 
to N, N sc being real, there is a t <-> (— t) symmetry in ( P7| ) and therefore 5^ = 2 J °° 
We now divide the time axis into four intervals : 

< t < At: The shortest time scale in our problem is At = h/AE. The 
contributions to this time interval are due to the differences between the exact 
and the semiclassical mean densities of states. This is an important observation, 
since it allows us to distinguish between the two sources of semiclassical error — the 
error that emerges from the mean densities and the error that originates from the 
fluctuating part (periodic orbits). Since we are interested only in the semiclassical 
error that results from the fluctuating part of the spectral density, we shall ignore 
this regime in the following. 

At < t < t crg : This is the non-universal regime ||29|| , in which periodic orbits are 



still sparse, and cannot be characterized statistically in a significant fashion. The 
"ergodic" time scale t eTg is purely classical and is independent of Ti. 

t cr g < t < tjj: In this time regime periodic orbits are already in the universal regime 
and are dense enough to justify a statistical approach to their proliferation and 
stability. The upper limit of this interval is the Heisenberg time £h — hd(E), which 
is the time that is needed to resolve the quantum (discrete) nature of a wavepacket 
with energy concentrated near E. The Heisenberg time is "quantum" in the sense 
that it is dependent of h: t H = 0{h l ~ d ). 

tri < t < oo: This is the interval of "long" orbits which is effectively truncated from the 
integration as a result of introducing a smoothing of the quantal and semiclassical 
counting functions, with a smoothing scale of the order of a mean level spacing. 

Dividing the integral (|27|) according to the above time intervals, we can rewrite Smooth : 

= ( f + r + r) 2 4 m E) - D U f, e) - 



h 

= ^short + + ^lorlg ■ (32) 

(2) 

As explained above, 5 lo ^ g can be ignored due to smoothing on the scale of a mean 

(2) 

level spacing. The integral <5 short is to be neglected for the following reason. The integral 
extends over a time interval which is finite and independent of h, and therefore it contains 
a fixed number of periodic orbits contributions. The semiclassical approximation 
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provides, for each individual contribution, the leading order in h, and therefore [[J(J 
we should expect: 



6. 



(2) 
short 



as ft 



0. 



(33) 



The purpose of this work is to check whether the semiclassical error is finite or divergent 
as h — ► 0, and to study if the rate of divergence depends on dimensionality. Equation 
(|33|) implies that 5^ OIt cannot affect 5^ in the semiclassical limit and we shall neglect 
it in the following. 

We thus remain with a lower bound for our measure: 



tf (2) S3 
"smooth 



5 m 

■ ■ . 



(34) 



which is going to be our object of interest from now on. 

The fact that in is extremely large on the classical scale renders the calculation of 
all the periodic orbits with periods less than tn an impossible task. However, sums over 
periodic orbits when the period is longer than t erg tend to meaningful limits, and hence, 
we would like to recast the expression for 8$ in the following way. Write 8$ as: 



5 (2) 



hi 



dt 



hi 



dt 



D(t) - D sc (t) 



D(t) 



D(t) - D sc {t) 



bit) 



hi 



dt 



D(t) ) x C{t) 



(35) 



(36) 



(37) 



/ envelope x correlation 



where the parametric dependence on E was omitted for brevity. The smoothing over t 
is explicitly indicated to emphasize that one may use a statistical interpretation of the 
terms of the integrand. This is so because in this domain, the density of periodic orbits 
is so large, that within a time interval of width h/AE there are exponentially many 
orbits whose contributions are averaged due to the finite resolution. 

We note now that we can use the following relation between the time spectrum and 
the spectral form factor K{r): 

2^ 



D(t) 



dt 



K(r) 



, , 9 dr (38) 

where r = t/t-R is the scaled time. The above form factor is smoothed according to the 
window function w. Hence: 



5 (2) S3 

smooth 



2tt 2 



dr 



K(t)C(t) 



T 



(39) 
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For generic chaotic systems we expect that K(t) agrees with the results of RMT in the 
universal regime r > r crg [[IT], [^, |2!| , and therefore 



K(t) ss gr for r crg < r < 1, (40) 

where <? = 1 for systems which violate time reversal symmetry, and g = 2 if time reversal 
symmetry is respected. This implies that the evaluation of <5 smooth reduces to 



Smooth- 2^2/ dT ^ _ • ( 41 ) 



r crst 



The dependence on ft in this expression comes from the lower integration limit which is 
proportional to h d ~ l as well as from the implicit dependence of the function C on h. 

Formula fl4l|) is our main theoretical result. However, we do not know how to 
evaluate the correlation function C(r) from first principles. The knowledge of the h 
corrections to each of the terms in the semiclassical time spectrum is not sufficient since 
the resulting series which ought to be summed is not absolutely convergent (see detailed 
discussion in section (|4j)). Therefore we have to recourse to a numerical analysis, which 
will be described in the next section. The numerical approach requires one further 
approximation, which is imposed by the fact that the number of periodic orbits with 
t < tn is prohibitively large. We had to limit the data base of periodic orbit to the 
domain t < t cpu with t erg <C t cpu <C t H - The time t cpu has no physical origin, and it 
represents only the limits of our computational resources. Using the available numerical 
data we were able to compute C(t) numerically for all t erg < t < t cpu and we then 
extrapolated it to the entire domain of interest. We consider this extrapolation procedure 
to be the main source of uncertainty. However, since the extrapolation is carried out 
in the universal regime, it should be valid if there are no other time scales between t eTg 
and £h- 



3. Numerical results 

We used the formalism and definitions presented above to check the accuracy of the 
semiclassical spectra of the 2D and 3D Sinai billiards. The most important ingredient 
in this numerical study is that we could apply the same analysis to the two systems, 
and by comparing them to give a reliable answer to the main question posed in this 
work, namely, how does the semiclassical accuracy depend on dimensionality. 

The classical dynamics in billiards depends trivially on the energy (velocity), and 
therefore the relevant parameter is the length rather than the period of the periodic 
orbits. Because of the same reason, the quantum wavenumbers k n = y/2mE n /h are 
the relevant variables in the quantum description. From now on we shall use the 
variables (l,k) instead of (t,E), and use "length spectra" rather than "time spectra". 
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Figure 4. Left: The quarter desymmetrized 2D Sinai billiard. Right: The fully (1/48) 
desymmetrized 3D Sinai billiard (bold lines). 



The semiclassical limit is obtained for k — > oo and 0(h) is equivalent to 0(k~ 1 ). Note 
also that for a billiard N(k) ~ Ak d where A is a proportional to the billiard's volume. 

The numerical work is based on the quantum spectra and on the classical periodic 
orbits which were calculated by Schanz and Smilansky |32| , |33| for the 2D billiard, and 
by Primack and Smilansky |5], [7| for the 3D billiard. The numerical methods and the 
checks performed to ensure that the quantum and the classical data bases are accurate, 
complete and immaculate are discussed in the papers cited above. 

We start with the 2D Sinai billiard, which is the free space between a square of 
edge L and an inscribed disc of radius R, with 2R < L. In our case we used L = 1 and 
R = 0.25 and considered the quarter desymmetrized billiard (see figure f|) with Dirichlet 
boundary conditions for the quantum calculations. The quantal data base consisted of 
the lowest 27645 eigenvalues in the range < k < 1320, with eigenstates which are 
either symmetric or antisymmetric with respect to reflection on the main diagonal. 
The classical data base consisted of the shortest 20273 periodic orbits (including time 
reversal, reflection symmetries and repetitions) in the length range < / < 5. For each 
orbit, the length, the stability determinant and the reflection phase were recorded. 

We begin the numerical analysis of the 2D Sinai billiard by numerically 
demonstrating the correctness of equation fl33|). That is, that for each individual periodic 
orbit, the semiclassical error indeed vanishes in the semiclassical limit. In figure |5] we 
plot \D — D sc \ for I = 0.5 as a function of k. This length corresponds to the shortest 
periodic orbit, that is, the one that runs along the edges that connect the circle with the 
outer square. For D sc we used the Gutzwiller trace formula. As is clearly seen from the 
figure, the quantal-semiclassical difference indeed vanishes (approximately as A; -1 ), in 
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k 



Figure 5. The absolute difference between the quantal and the semiclassical 
(Gutzwiller) length spectra for the 2D Sinai billiard at I = 0.5. This length corresponds 
to the shortest unstable periodic orbit. The average log-log slope is about —1.1, 
indicating approximately decay. The data were averaged with a Gaussian window. 



accordance with fl33|). We emphasize again, that this behavior does not imply that 5^ 
vanishes in the semiclassical limit, since the number of terms depends on k. It implies 

(2) 

only that o^ ort vanishes in the limit, since it consists of a fixed and finite number of 
periodic orbit contributions. We should also comment that (non-generic) penumbra 
corrections to individual grazing orbits introduce errors which are of order k~"< with 
< 7 < 1 j34], [13| . However, since the definition of "grazing" is in itself k dependent, 

(2) 

one can safely neglect penumbra corrections in estimating the large k behavior of £ s ho rt . 

We now turn to the main body of the analysis, which is the evaluation of 5$ for 
the 2D Sinai billiard. Based on the available data sets, we plot in figure ^| the function 
C(l; k) in the interval 2.5 < / < 5 for various values of k. One can observe, that as a 
function of / the functions C(l; k) fluctuate in the interval for which numerical data were 
available, without exhibiting any systematic mean trend to increase or to decrease. We 
therefore approximate C(l; k) by 



C(l; k) « const • f(k) = C avg (k). 



(42) 
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According to the discussion in section |2| we extrapolate this formula in I up to the 
Heisenberg length L H = 2nd(k) and using fl4T|) we obtain: 

HZ = hi( WW = C avg (k) 0(ln fc). (43) 

The last equality is due to L H = C , (fc a! ~ 1 ). To evaluate C avg (fc) we averaged C(l;k) 
over the interval L erg = 3.5 < / < 5 = L cpn and the results are shown in figure [7[ We 
choose L erg = 3.5 because the density of periodic orbits is already large for this length 
(see figure ^) and we expect universal behavior of the periodic orbits. (For the Sinai 
billiard described by flow the approach to the invariant measure is algebraic rather than 
exponential ]35], p6| , and thus one cannot have a well-defined L erg . An any rate, the 
specific choice of L erg did not affect the results in any appreciable way.) Inspecting 
C avg (A;), it is difficult to arrive at firm conclusions, since it seems to fluctuate around 
a constant value up to k « 900 and then to decline. If we approximate C avg (A;) by a 
constant, we get a "pessimistic" value of <5^: 

*£S,(*) = k ) = h ) "pessimistic" (44) 
while if we assume that C avg (fc) decays as a power-law, C avg (/c) = k~^,f3 > 0, then 

CShW = 0(k~P Ink) — > "optimistic" . (45) 
Collecting the two bounds we get: 

0{k-f> In k) < «52SL(fc) < 0(ln fc) • (46) 

Our estimates for the 2D Sinai billiard can be summarized by saying that the 
semiclassical error diverges no worse than logarithmically (meaning, very mildly) . It may 
well happen that the semiclassical error is constant or even vanishes in the semiclassical 
limit. To reach a conclusive answer one should invest exponentially larger amount of 
numerical work. 

There are a few comments in order here. Firstly, the quarter desymmetrization of the 
2D Sinai billiard does not exhaust its symmetry group, and in fact, a reflection symmetry 
around the diagonal of the square remains. This means, that the spectrum of the quarter 
2D Sinai billiard is composed of two independent spectra, which differ by their parity 
with respect to the diagonal. If we assume that the semiclassical deviations of the 
two spectra are not correlated, the above measure is the sum of the two independent 
measures. It is plausible to assume also that both spectra have roughly the same 
semiclassical deviation, and thus Smooth i s twice the semiclassical deviation of each 
of the spectra. Secondly, we recall that the 2D Sinai billiard contains "bouncing-ball" 



families of neutrally stable periodic orbits [37, p8l B2J. We have subtracted their leading- 



order contribution from D such that it includes (to leading order) only contributions 
from generic, isolated and unstable periodic orbits. This is done since we would like to 
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Figure 6. The functions C(l;k) for quarter 2D Sinai billiard L = 1,R= 0.25 with 
Dirichlet boundary conditions. The window w(k' — k) was taken to be a Gaussian with 
standard deviation a = 60. We averaged C(l; k) over Z-intervals of « 0.2 in accordance 
with ( p6| ) to avoid sharp peaks due to small denominators. The averaging, however, 
is fine enough not to wash out all of the features of C(l; k). The vertical bars indicate 
the locations of primitive periodic orbits, and the daggers indicate the locations of the 
bouncing-ball families. 



deduce from the 2D Sinai billiard on the 2D generic case in which the bouncing-balls 
are not present. (In the Sinai billiard, which is concave, there are also diffraction effects 



34] , [0|, but we did not treat them here.) Thirdly, the variant of (|38|) for billiards reads: 

^ ( °r), d/= Sl de (47) 

when £ = I/L-r. In figure |] we demonstrate the compliance of the form factor with 
RMT GOE using the integrated version of the above relation, and taking into account 
the presence of two independent spectra. Fourthly, it is interesting to know the actual 

(2) 2D 

numerical values of Smooth W f° r the k values that we considered. We carried out the 
computation, and the results are presented in figure |9|. It is interesting to observe, that 

(2) 2D 

for the entire range we have Smooth (^) ~ 0-1 ^ 1> which is very encouraging from an 
"engineering" point of view. 

We now turn to the analysis of the 3D Sinai billiard. The billiard is the free space 
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Figure 7. Averaging in I of C(l; k) for 2D Sinai billiard as a function of k. 



between a cube of edge L and an inscribed sphere of radius R, where 2R < L (see figure 
^). We used L = 1 and R = 0.2 and desymmetrized the billiard to the fundamental 
element (1/48 of the original one). We calculated the lowest 6697 quantum levels in the 
interval < k < 281.1 and the shortest 586965 periodic orbits with length < / < 5 
(the number includes repetitions and time-reversal conjugates). 

To treat the 3D Sinai billiard we need to somewhat modify the formalism which was 
presented in the previous section. This is due to the fact that in 3D the contributions 
of the various non-generic bouncing-ball manifolds overwhelm the spectrum || [7|. 
Since our goal is to give an indication of the semiclassical error in generic systems, 
it is imperative to get rid of this very strong non-generic effect. The bouncing ball 
amplitudes are 0{k^ s ~ 1 ^ 2 ) where s is the dimensionality of bouncing-ball manifold 
in configuration space. In 3D s = 3 typically, which completely overwhelms the 
contributions from isolated periodic orbits whose amplitude is O(k ). Even the 
diffraction corrections to the bouncing-ball amplitudes in 3D increase as 0{k 1 ) with 
7 > 0. In contrast with the 2D case, however, it is difficult to subtract the bouncing- 
ball contributions analytically for two main reasons. First, there are always infinitely 
many bouncing-ball primitive families in the 3D case, while in 2D there is only a 
finite number. Indeed, in any finite length interval there exists only a finite number 
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Figure 8. Verification of equation for the quarter 2D Sinai billiard. We plot 
/(£) = d£' K(£')/{;' 2 and compare the quantum data with RMT. The minimal £ 
corresponds to L crg = 3.5. The integration is done for smoothing, and we fix the upper 
limit to avoid biases due to non-universal regime. Note the logarithmic scale. 



of bouncing ball lengths, however their number for billiards in 3D exceeds by far the 
corresponding number for billiards in 2D. Second, the semiclassical amplitudes of the 
bouncing-balls are proportional to their volume in configuration space, and it is difficult 
in general to calculate it analytically. Thus, to overcome these difficulties we had to 
device a special method to cleanse the spectrum from the effect of the bouncing-balls 
and from the leading diffractive corrections. This method relies on the sensitivity of the 
eigenvalues to the boundary condition and it is described in detail in [7|]. We shall 
briefly describe the essence of the method. 

The most general ( "mixed" ) boundary conditions under which the quantum billiard 



problem is self-adjoint can be written as |21], |39 



k, cos a tp(r ) + sin a df{ip(r ) = , r G (boundary of the billiard) (48) 

where n is the normal pointing outside of the billiard, the angle a interpolates smoothly 
between Dirichlet (a = 0) and Neumann (a = 7r/2) cases and k has the dimension of a 
wavenumber. Note that a and k can be different for different parts of the boundary. The 
spectrum depends parametrically on the boundary parameters, {k n (a, k)}, and thus we 
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Figure 9. The numerical values of #gJ ooth for the quarter 2D Sinai billiard. We 

(2) 

included also the contribution <r, ; of the non-universal regime. The contributions 

(2) (2) 

from the time interval t erg < t < t cpu are contained in $m,cpu, and 8^ ext is the 
extrapolated value for t cpu < t < in (refer to ([32]) and to the end of section ||) . 



can define the counting function N(k; a, k). In |^9|, 0] we discussed the 2D and 3D Sinai 
billiards, and showed that if we choose Dirichlet boundary conditions on all the billiard 
boundaries excluding the circle (sphere), and set k ^ on the circle (sphere) then the 
quantity 

a=0 n 

is to a large extent free of the effects of the bouncing-balls. In the above k n = k n (a = 0) 
which are the Dirichlet-everywhere eigenvalues. The corresponding semiclassical trace 
formula reads [[39], 0: 

d sc {k] k) = (smooth part) + AjBj cos(kLj — Vjir/2) (50) 

po 

where 



d(k; k) = 



dN(k; a, k) 
da 



E 



r>=0 



dk n (a, k) 
da 
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In the above rij is the number of collisions with the circle (sphere) of the j'th periodic 
orbit, and {e[ j) } are the angles of incidence on the boundary with respect to the normal. 
We note that ( fj9|) is a weighted density of states where the standard 1 weights of the 5 
functions are replaced by v n , and (|50| ) is a weighted trace formula where the standard 
amplitudes Aj are replaced by AjBj. One can show, that Bj « v n for long enough 
(ergodic) orbits. 

We shall use d for our purposes as follows. Let us consider the weighted counting 
function: 

N(k) = t dk' d(k') = Yv n e(k - k n ). (52) 

The function AT is a staircase with stairs of variable height v n . As was explained above, 
its advantage over A" is that it is semiclassically free of bouncing balls effects (to leading 
order) and corresponds only to the generic periodic orbits p9| . Similarly, we construct 
from d sc the function A" sc . Having defined N, N sc , we proceed in analogy to the Dirichlet 
case. We form from N, N sc the functions N, N sc , respectively, by multiplying with a 
window function w(k' — k) and then construct the measure 

5 (2) 

as in (E5J) . The only 

difference is that the normalization of w must be modified to account for the "velocities" 
v n such as: 

d-\k)J2<Hkn-k)\ 2 = l. (53) 
n 

The above considerations are meaningful provided that the "velocities" v n are narrowly 
distributed around a well-defined mean v(k) and we consider a small enough fc-interval, 
such that v(k) does not change appreciably within this interval. Otherwise, is 
greatly affected by the fluctuations of v n (which is undesired) and the meaning of the 
normalization is questionable. 

To demonstrate the utility of the above construction using the mixed boundary 
conditions, we return to the 2D case. We set k = 1007T, and note that the spectrum at 
our disposal for the mixed case was confined only to the interval < k < 600. First 
of all, we want to examine the width of the distribution of the t> n 's. In figure |10| we 



plot the ratio of the standard deviation of v n to the mean, averaged over the fc-axis 
using a Gaussian window. We use the same window also in the calculations below. The 
observation is that the distribution of v n is moderately narrow and the width decreases 
algebraically as k increases. This justifies the use of the mixed boundary conditions as 
was discussed above. One also needs to check the validity of (pETD, and indeed we found 
compliance with GOE also for the mixed case (results not shown). We next compare 
the functions C(l; k) for both the Dirichlet and the mixed boundary conditions. It turns 
out, that also in the mixed case the functions C(l; k) (not shown) fluctuate in / with no 
special tendency. The averages C ecvg {k) for the Dirichlet and mixed cases are compared 
in figure |TT|. The values in the mixed case are systematically smaller than in the Dirichlet 
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case which is explained by the efficient filtering of tangent and close to tangent orbits 



that are vulnerable to large diffraction corrections [33, 13]. However, from k = 250 on 
the two graphs show the same trends, and the values of C avg in both cases are of the 

(2) 

same magnitude. Thus, the qualitative behavior of <5 smooth is shown to be equivalent in 

(2) 

the Dirichlet and mixed cases, which gives us confidence in using Smooth together with 
the mixed boundary conditions procedure. 

We finally applied the mixed boundary conditions procedure to compute £ smoo th for 
the desymmetrized 3D Sinai with L = 1, R = 0.2 and set k = 100. We first verified that 



also in the 3D case the velocities v n have narrow distribution — see figure [TO]. Next, 
we examined equation (|47| ) using quantal data, and discovered that there are deviations 
form GOE (figure |12|). We have yet no satisfactory explanation of these deviations, 
but we suspect that they are caused because the ergodic limit is not yet reached for 
the length regime under consideration due to the effects of the infinite horizon which 
are more acute in 3D. Nevertheless, from observing the figure as well as suggested by 
semiclassical arguments, it is plausible to assume that K(£) oc £ for small £. Hence, this 
deviation should not have any qualitative effect on 5^ according to (]4T|). Similarly to the 
2D case, the behavior of the function C(l; k) is fluctuative in I, with no special tendency 



(figure [13]). If we average C(l; k) over the universal interval L erg = 2.5 < / < L 



cpu 



we obtain C avg (/c) which is shown in figure [14j. The averages C avg (k) are fluctuating 



with a mild decease in k, and therefore we can conclude that 

0(k^lnk)<5^ h <0(lnk) (54) 

where the "optimistic" measure (leftmost term) corresponds to C avg ( k) = C(/c _/3 ),/3 > 0, 
and the "pessimistic" one (rightmost term) is due to C avg (/c) = const. In other words, 
the error estimates fl4j| HD for the 2D and the 3D cases, respectively, are the same, and 



in sharp contrast to the "traditional" error estimate which predicts that the errors should 
be different by a factor 0(h^ 1 ). On the basis of our numerical data, and in spite of 
the uncertainties which were clearly delineated, we can safely exclude the "traditional" 
error estimate. 



4. Discussion 



Our main finding was that the upper bound on the semiclassical error is a logarithmic 
divergence, which is independent of the dimension (equations (f46|), (|54])). In this respect, 
there few points which deserve discussion. 

(2) 

To begin, we shall try to evaluate <5 smooth using the explicit expressions for the 
leading corrections to the semiclassical counting function of 2D generic billiard system, 
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Figure 10. Calculation of Q = y/yu%) — { v n) 2 /\{v n )\ for quarter 2D Sinai billiard 
(up) and for the desymmetrized 3D Sinai billiard (down). 



as derived by Alonso and Gaspard JTT[: 



A; 



N(k) = N(k) + J2 — sin 



% + 0(l/k 2 ) 
k 



(55) 
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k 

Figure 11. Comparison of C &vg (k) for Dirichlet and mixed boundary conditions for 
the quarter 2D Sinai billiard. We used a Gaussian window with a = 40. 



where A,- are the standard semiclassical amplitudes (see fl30D), Lj are the lengths of 
periodic orbits and Qj are the fc-independent amplitudes of the 1/k corrections. The 
Q/s are given in . We ignored in the above equation the case of odd Maslov indices. If 
we calculate from N(k) the corresponding length spectrum D(l; k) using a (normalized) 
Gaussian window w{k' — k) = (l/\Ar(7 2 ) exp[— (k' — k) 2 /(2a 2 )], we obtain: 



eifc( i_ L .)_ i ^ e _ a _ Li) 2^ _ Qik{l+Lj)+i S± Q _ {l+L] f S l 



.(56) 



In the above we regarded the phase e^ J//fe as slowly varying. The results of Alonso and 
Gaspard |TTJ suggest that the Qj are approximately proportional to the length of the 
corresponding periodic orbits: 

Qj w QLj. (57) 

We can therefore well-approximate D as: 

D(l; k) » ^e~^ k Y,T~ [•••]= e" tQi/fc 4c- GTF , (58) 
Z V 71 j L j 

where D sc _gtf is the length spectrum which corresponds to the semiclassical Gutzwiller 
trace formula for the counting function (without 1/k corrections). We are now in a 
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Figure 12. Check of equation (E^) for the desymmetrized 3D Sinai billiard. The 
minimal £ corresponds to i crg = 2.5. The function is defined as in figure ^. Note 
the logarithmic scale. 



position to evaluate the semiclassical error, indeed: 



CootM = 2 1 " dl D(l; k) - Ac-gtf(Z; k) 



-^mii 



b{i-k) 



(59) 



If we now use equation (f47|) and K(l) « gl/L^ (which is valid for / < Lh for chaotic 
systems), we get: 



^smooth (k) 



2g [Lu dl . 2 (Ql\ 2g rQLn/m sin 2( f ) 



ir 2 Jr ■ I 



sin 



, , dt- 

2k J 7l 2 jQL IIli J(2k) t 



For fc^oowe have that 

r QL min /(2k) s i n 2 (t) rQL min /(2k) 



dt- 



t 



dt-t = 0{l/k 2 



which is negligible, hence we can replace the lower limit in (|60|) with 0: 



^smooth (ty 



2g 



7T 2 JO 



~2k~ sin [t) 



t 



(60) 



(61) 



(62) 
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Figure 13. The functions C(l; k) for desymmetrized 3D Sinai billiard L = 1, R = 0.2 
with mixed boundary conditions. We took a Gaussian window with a = 20, and 
smoothed over /-intervals of « 0.3. The upper vertical bars indicate the locations of 
primitive periodic orbits. 



This is the required expression. The dimensionality enter in £ smoo th(^) only through the 
power of k in Lh- 

Let us apply equation ( |52j ) to 2D and 3D cases. For 2D we have that Lh = Ak in 
leading order, where A is the billiard's area, thus, 

«££-(*) - | J™' = oonst = 0(k°) (63) 

which means, that the semiclassical error in 2D billiards is of the order of the mean 
spacing, and therefore the semiclassical trace formula is (marginally) accurate and 
meaningful. This is compatible with our numerical findings within the limitations of 
the numerical fluctuations. 

For 3D, the coefficients Qj were not obtained explicitly, but we shall assume that 
they are still proportional to Lj (equation (|57|)) and thus that fl62|) holds. For 3D 
billiards Lh = (V/7r)k 2 to leading order, where V is the billiard's volume. Thus the 
upper limit in ( |62| ) is QVk/(2ir) which is large in the semiclassical limit. In that case, 
we can replace sin 2 (t) with its mean value 1/2 and the integrand becomes essentially 
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Figure 14. Averaging in I of C(l;k) for 3D Sinai billiard as a function of k. The 



averaging was performed in the interval L 



erg 



2.5 < I < 5 = L. 



cpu- 



1/t which results in: 

*2£La(*) = 0(lnfc). 



(64) 



That is, in contrast to the 2D case, the semiclassical error diverges logarithmically and 
the semiclassical trace formula becomes meaningless as far as the prediction of individual 
levels is concerned. This is compatible with our numerical results within the numerical 
dispersion. However, it relies heavily on the assumption that Qj ~ QLj, for which we 
can offer no justification. We note in passing, that the logarithmic divergence persists 
also for d > 3. 

Another interesting point relates to integrable systems. It can happen, that for an 
integrable system it is either difficult or impossible to express the Hamiltonian as an 
explicit function of the action variables. In that case, we cannot assign to the levels 
other quantum numbers than their ordinal number, and the semiclassical error can be 
estimated using S^ 2 \ However, since for integrable systems K(t) = 1, we get that: 

,( 2) ,int _ J_ f 1 , C{t) 

^smooth ~ J T aT ^2^- 

Therefore, for deviations which are comparable to the chaotic cases, C(r) = 0(1), we 
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get Smooth = 0{h l ~ d ) which is much larger than for the chaotic case and diverges for 
d> 1. 

The formula (|T) for the semiclassical error contains semiclassical information in 
two respects. Obviously, C(r), which is the difference between the quantal and the 
semiclassical length spectra contains semiclassical information. But also, the fact that 
the lower limit of the integral in ( f4T|) is finite is a consequence of semiclassical analysis. 
If this lower limit is replaced by 0, the integral diverges for finite values of h, which is 
meaningless. Therefore, the fact that the integral has a lower cutoff, or rather, that D is 
exactly below the shortest period, is a crucial semiclassical ingredient in our analysis. 

Finally, we consider the case in which the semiclassical error is estimated with no 
periodic orbits taken into account. That is, we want to calculate (\N(E) — N(E)\ 2 )e 
which is the number variance E 2 (x) for the large argument x = AE d(E) ^> 1. This 
implies C(r) = 1, and using (f4~l"D we get that 5^ ooth = g/(2ir 2 ) ln(t H Aerg), which in the 
semiclassical limit becomes g/(2ir 2 ) ln(tn) = 0(lnh). This result is fully consistent and 
compatible with previous results for the asymptotic (saturation) value of the number 
variance S 2 (see for instance [40], 41|). It implies also that the pessimistic error 



bound (f44|) is of the same magnitude as if periodic orbits were not taken into account 
at all. (Periodic orbits improve, however, quantitatively, since in all cases we obtained 
Cavg < 1-) Thus, if we assume that periodic orbit contributions do not make N sc worse 
than N, then the pessimistic error bound 0(lnh) is the maximal one in any dimension 
d. This excludes, in particular, algebraic semiclassical errors, and thus refutes the 
traditional estimate 0(h 2 ~ d ). 
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